from typing import Any, Dict, Optional, Tuple, cast

import numpy as np
from icecream import ic

from .fitmodel import (
    MeshEnum,
    MidTailDcttRet,
    MidTailDtctCfg,
    MidTailDtctData,
    NmlzIn,
    WtdStdDevIn,
)
from .fitmodel import myDftVals as dv


##0000000000000000000000000000000000000000000000000000000000000000000000000 拟合工具箱
def normalize_by_sum(nIn: NmlzIn) -> np.ndarray:
    """
    输入 arr1D, 除以 sum 作归一化; 再除以 area, 从 分 bin dose -> density 分布

    Parameters:
    - mtype: mesh type; box, cylinder
    - arr1D: Input array (dose values)
    - binwidth: Bin width for normalization

    Returns:
    - Normalized array (density distribution)
    """
    total_sum = np.sum(nIn.vars_tgt)
    if total_sum == 0:
        return np.zeros_like(nIn.vars_tgt)

    if nIn.mtype == MeshEnum.Box:
        return nIn.vars_tgt / total_sum / nIn.bin_width

    if nIn.mtype == MeshEnum.Cyliner:
        rs = np.arange(len(nIn.vars_tgt), dtype=float)  # 创建一个索引数组
        rs *= nIn.bin_width
        areas = np.pi * (1 + 2 * rs)  ## 圆环面积
        areas[0] = np.pi * np.pow(nIn.bin_width, 2.0)  ## 中心圆盘面积
        return nIn.vars_tgt / areas / total_sum

    raise ValueError(f"wrong mesh type:{nIn.mtype}")


def weighted_stddev(wIn: WtdStdDevIn) -> float:
    """
    Calculate the weighted standard deviation of x with weights y.

    Parameters:
    - qty: Input values
    - wts: Weights (typically a density or probability distribution)
    - eps: Small value to prevent division by zero

    Returns:
    - Weighted standard deviation
    """
    wts = np.maximum(wIn.wts, wIn.eps)
    wts_tot = np.sum(wts)
    if wts_tot == 0:
        return 0.0
    iwts = wts / wts_tot
    mu = np.average(wIn.qty, weights=iwts)
    var = np.average(np.pow(wIn.qty - mu, 2.0), weights=iwts)
    return np.sqrt(var)


def mid_tail_weights(mIn: MidTailDcttRet) -> np.ndarray:
    """
    Exponential weighting function to emphasize regions around the center value.

    Parameters:
        x: Input coordinates
        center: Center of the weighting region
        width: Width of the weighting region
        lam_m: Weighting strength

    Returns:
        Weight array
    """
    crd_abs = np.abs(mIn.crds)
    iExpo = -0.5 * np.pow((crd_abs - mIn.center) / max(mIn.width, dv.epsilon), 2.0)
    return 1.0 + mIn.lam_m * np.exp(iExpo)


def mid_tail_detect_v2(
    mIn: MidTailDtctData,
    cfg: Optional[MidTailDtctCfg] = None,
) -> MidTailDcttRet:
    """
    Detect regions where the predicted values are significantly below the actual values,
    typically in the mid-tail regions of the distribution.

    Parameters:
        crds: Input coordinates
        vals: Actual values
        vals_pred0: Initial prediction values
        config: MidTailConfig object with detection parameters

    Returns:
        Dictionary containing center, width, and lam_m parameters
    """
    # Use default config if not provided
    if cfg is None:
        cfg = MidTailDtctCfg()

    # r_min = config.r_min
    # r_max = cfg.r_max
    # smooth_window = cfg.smooth_window
    # width_clip = cfg.width_clip
    # lam_m_clip = cfg.lam_m_clip
    # rel_err = cfg.rel_err
    # abs_err = cfg.abs_err
    vals_diff = mIn.vals_p0 - mIn.vals_tgt
    eps = dv.epsilon
    threshold = np.maximum(cfg.abs_err, cfg.rel_err * np.maximum(mIn.vals_tgt, eps))

    # Find regions where y_pred < y by more than threshold
    mask1 = vals_diff < -threshold
    mask2 = (np.abs(mIn.crds) >= cfg.r_min) & (np.abs(mIn.crds) <= cfg.r_max)
    ix_mask = mask1 & mask2
    if not np.any(ix_mask):
        return MidTailDcttRet(
            crds=mIn.crds,
            center=dv.mtl_center,
            width=5.0,
            lam_m=1.0,
            left=None,
            right=None,
        )

    crd_s = mIn.crds[ix_mask]  ## coords selected
    ## relative value difference
    rvd_s = -(vals_diff[ix_mask]) / np.maximum(mIn.vals_tgt[ix_mask], eps)

    # Sort by x values
    ix_sorted = np.argsort(crd_s)
    crd_st = crd_s[ix_sorted]
    rvd_st = rvd_s[ix_sorted]

    # Find longest continuous region where y_pred < y
    if len(crd_st) < 3:
        return MidTailDcttRet(
            crds=mIn.crds,
            center=cast(float, np.mean(crd_st)),
            width=2.0,
            lam_m=8.0,
            left=crd_st.min(),
            right=crd_st.max(),
        )

    gap_threshold = dv.gap_threshold  # 1.5 * dx, where dx is 0.1

    # Find start and end points of continuous regions using numpy operations
    # Calculate differences between consecutive x values
    crd_diffs = np.diff(crd_st)
    # Find gap locations where difference > threshold
    gap_indices = np.where(crd_diffs > gap_threshold)[0]
    # Add boundaries: start at 0, end at len(x_sorted)
    breakpoints = np.concatenate(([0], gap_indices + 1, [len(crd_st)]))

    ## Find the longest continuous region; segment
    seg_lengths = np.diff(breakpoints)
    k_max = int(np.argmax(seg_lengths))
    left, right = breakpoints[k_max], breakpoints[k_max + 1]
    crd_seg, rvd_seg = crd_st[left:right], rvd_st[left:right]

    if len(crd_seg) < 3:
        return MidTailDcttRet(
            crds=mIn.crds,
            center=cast(float, np.mean(crd_seg)),
            width=2.0,
            lam_m=8.0,
            left=crd_seg.min(),
            right=crd_seg.max(),
        )

    ## Smooth the relative values difference
    if cfg.smooth_window > 1 and cfg.smooth_window < len(rvd_seg):
        kernel = np.ones(cfg.smooth_window) / cfg.smooth_window
        rvd_sm = np.convolve(rvd_seg, kernel, mode="same")
    else:
        rvd_sm = rvd_seg

    # Find peak of the relative error
    ix_peak = np.argmax(rvd_sm)
    crd_peak = crd_seg[ix_peak]
    rvd_peak = rvd_sm[ix_peak]
    print(f"peak_val: {rvd_peak:.3f}")

    # Calculate width as half maximum width
    rvd_half = 0.1 * rvd_peak

    # Find leftmost index where value drops below half using numpy
    rvd_to_left = rvd_sm[: ix_peak + 1][::-1]  # Reverse order from peak to left
    left_indices = np.where(rvd_to_left < rvd_half)[0]
    iL = ix_peak - left_indices[0] if len(left_indices) > 0 else 0

    # Find rightmost index where value drops below half using numpy
    rvd_to_right = rvd_sm[ix_peak:]
    right_indices = np.where(rvd_to_right < rvd_half)[0]
    iR = ix_peak + right_indices[0] if len(right_indices) > 0 else len(crd_seg) - 1

    x_left, x_right = crd_seg[iL], crd_seg[iR]
    width = 0.5 * abs(x_right - x_left)
    width = np.clip(width, cfg.width_clip[0], cfg.width_clip[1])

    # Calculate lambda_m parameter
    lam_m = 8.0 + 30.0 * (1 - np.exp(-3.0 * rvd_peak))
    lam_m = np.clip(lam_m, cfg.lam_m_clip[0], cfg.lam_m_clip[1])

    return MidTailDcttRet(
        crds=mIn.crds,
        center=crd_peak,
        width=width,
        lam_m=lam_m,
        left=x_left,
        right=x_right,
    )
